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Abstract: The DFLU numerical flux was introduced in order to solve hyperbolic scalar conservation laws 
with a flux function discontinuous in space. We show how this flux can be used to solve systems of 
conservation laws. The obtained numerical flux is very close to a Godunov flux. As an example we 
consider a system modeling polymer flooding in oil reservoir engineering. 
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Application du flux DFLU aux systemes de lois de conservation 



Resume : Le flux numerique DFLU a ete introduit afin de resoudre des lois de conservation scalaires 
hyperboUque avec des fonctions de flux discontinues en espace. Nous montrons comment ce flux peut 
etre utilise pour resoudre des systemes de lois de conservation. On obtient ainsi un flux numerique tres 
proche du flux de Godunov. Comme exemple on considere un systeme modelisant I'injection de polymere 
en ingenierie de reservoir petrolier. 

Mots-cles : Volumes finis, differences finies, solveurs de Riemann, systemes de lois de conservation, 
ecoulements en milieu poreux, injection de polymeres. 
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1 Introduction 

The main difficulty in the numerical solution of systems of conservation laws is the complexity of con- 
structing the Riemann solvers. One way to overcome this difficulty is to consider centered schemes as in 
IIT5I ITSl I22I I23I O. However, in general these schemes are more diffusive than Godunov type methods 
based on exact or approximate Riemann solvers when this alternative is available. Therefore in this paper 
we will consider Godunov type methods. Most often the numerical solution requires the calculation of 
eigenvalues or eigenvectors of the Jacobian matrix of the system. This is even more complicated when the 
system is non-strictly hyperbolic, i.e. eigenvectors are not linearly independent. In this paper we present a 
new approach which do not require such eigenvalue and eigenvector calculations. 
Let us consider a system of conservation laws in conservative form 

Ut + (F(U)), = 0, \5^{u\--- ,u-'), F= (/I, ••■ ,/■'). 

A conservative finite volume method reads 



TTn+l Tjn F" , - F" , 
At h 



= 



where F"_^]^/2 ^ numerical flux calculated using an exact or approximate Riemann solver. In a first order 
scheme this numerical flux is calculated using the left and right values U" and U"_^_]^. If we solve the 
equation field by field the j-th equation reads 



At h 
where the j-th numerical flux is a function of U" and 



= 



77'?:'"' 7717/ l,""- J,''^ 1-^ J.n \ ■ 1 7 

^/+i/2 = ^K 3^l,---,J. 
This flux function can be calculated by solving the scalar Riemann problem for t > tn'. 

Ui + if'^"{u^,x)):,^0, (1) 
U^{x,tn) = U-'" if a; < U^X,tn) = uf^^ if X > Xi+i/2, 

where the flux function , discontinuous at the point x = x 1^112, is defined by 



(2) 



(L and R refer to left and right of the point Xi^i/2)- 

Scalar conservation laws like equation ([T]l with a flux function discontinuous in space have been the 
object of many studies ||2l[l2lLl4 8, 10, 13, 24, 25 6, 20 2, 16|. In particular, in |2| a Godunov type finite 
volume scheme was proposed and convergence to a proper entropy condition was proved, provided that the 
left and right flux functions have exactly one local maximum and the same end points (the case where the 
flux functions has exactly one local minimum can be treated by symetry). At the discontinuity the interface 
flux, that we call the DFLU flux, is given by the formula 



{ul,ur) = min|/L(min{uL,6'L}),/i?,(max{ufl,6'i?})|, (3) 



if / denotes the scalar flux function and 9l =argmax(/i), 9ji =argmax(/fl). When /l = /_b this formula 
is equivalent to the Godunov flux so formula ^ can be seen as an extension of the Godunov flux to the 
case of a flux function discontinuous in space. In the case of systems formula (O can be applied to the 
fluxes and /|^". 

To illustrate the method we consider the system of conservation laws arising for polymer flooding in 
reservoir simulation which is described in section |2l This system, or similar systems of equations, is non- 
strictly hyperbolic and is studied in several papers 11211 [121 [m i9il. For example in [12J the authors solve 
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Riemann problems associated to this system when gravity is neglected and therefore the fractional flow 
function is an increasing function of the unknown. In this case, the eigenvalues of the corresponding Jaco- 
bian matrix are positive and hence it is less difficult to construct Godunov type schemes which turn out to 
be upwind schemes. When the above model with gravity effects is considered, then the flux function is not 
necessarily monotone and hence the eigenvalues can change sign. This makes the construction of Godunov 
type schemes more difficult as it involves exact solutions of Riemann problems with a nonmonotonous 
fractional flow function. Therefore in section[3]we solve the Riemann problems in the general case when 
gravity terms are taken into account so the flux function is not anymore monotone. This will allow to 
compare our method with that using an exact Riemann solver In section |4] we consider Godunov type 
finite volume schemes. We present the DFLU scheme for the system of polymer flooding and compare 
it to the Godunov scheme whose flux is given by the exact solution of the Riemann problem. We also 
present several other possible numerical fluxes, centered like Lax-Friedrichs or Force, or upstream like the 
upstream mobility flux commonly used in reservoir engineering ||4]|5][l6l. Finally in section|5]we compare 
numerically the DFLU method with these fluxes. 



2 A system of conservation laws modeling polymer flooding 

A polymer flooding model for enhanced oil recovery in petroleum engineering was introduced in lfT9ll as 
the following 2x2 system of conservation laws 

st + fis,c).^ = 
(sc + a(c))t + (c/(s,c)), - 

where t > and x G M., {s,c) € I x I with / = [0, 1]. s — s{x, t) denotes the saturation of the wetting 
phase, so 1 — s is the saturation of the oil phase, c = c{x, t) denotes the concentration of the polymer in 
the wetting phase which we have normalized. Here the porosity was set to 1 to simplify notations. The flux 
function / is the Darcy velocity of the wetting phase ipi and is determined by the relative permeabilities 
and the mobilities of the wetting and oil phases, and by the influence of gravity: 

f{s, c)^ipi= '^N^f'w + - 52)A2(s, c)]. (5) 

The quantities Xi,£ ~ 1, 2 are the mobilities of the two phases, with £ = 1 referring to the wetting phase 
and £ — 2 referring to the oil phase: 

Kkri{s) 

^e{s, c) = T^' * ^ 1' 2' 

Mi(c) 

where K is the absolute permeability, and kre and p^ are respectively the relative permeability and the 
viscosity of the phase £. kri is an increasing function of s such that kri (0) = while kr2 is a decreasing 
function of s such that fcr2(l) — 0. Therefore A^, £ = 1,2 satisfy 

Ai = Ai (s, c)is an increasing functions of s, Ai(0,c) = Vc £ [0,1], 
A2 = A2(s, c) is a decreasing functions of s, A2(l, c) = Vc G [0, 1]. 

The idea of polymer flooding is to dissolve a polymer in the injected water in order to increase the viscosity 
of the injected wetting phase. Thus the injected wetting phase will not be able to bypass oil so one obtains 
a better displacement of the oil by the injected phase. Therefore pi{c) is increasing with c while p2 will 
be taken as a constant assuming there is no chemical reaction between the polymer and the oil. Therefore 
/ will decrease with respect to c. The function a = a(c) models the adsorption of the polymer by the rock 
and is increasing with c. 

Lp is the total Darcy velocity, that is the sum of the Darcy velocities of the two phases Lpi and (f2'- 
ip^ipi+ip2, fi^ [y+ (gi -g2)A2], '/?2 ^ , ^\ VP+{g2 -ffi)Ai]. 

Ai -\- A2 Ai + A2 
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<y9 is a constant in space since we assume that the flow is incompressible. The gravity constants gi , g2 of 
the phases are proportional to their density. 
To equation (|4]i we add the initial condition 

(s(a;,0),c(x,0)) = {so{x) , co{x)) . (7) 

Since the case when / is monotone was already studied in lfT2l[TTl . we concentrate on the nonmonotone 
case which is more complicated and corresponds to taking into account gravity. Therefore we assume that 
(y5 = so for the nonlinearities of the system (|4]l. We will assume also that phase 1 is heavier than phase 2 
(<?i > 52) so we can assume the following properties: 

(i) /(s, c) > 0, /(O, c) = /(I, c) = for all c e /. 

(ii) The function s f{s, c) has exactly one global maximum in / with 6 =argmax(/). 

(iii) fc{s, c) < V s e (0, 1) and for aU c e / 

(iv) The adsorbtion term a ~ a{c) satisfies 

a(0) = 0, /i(c) = ^(c) > 0, 44(c) < for aU c G /. 
etc dc' 

Typical shapes of functions / and a are shown in Fig. [T] We expand the derivatives in equations (|4]l and we 




plug the resulting first equation into the second one. Then we obtain the system in nonconservative form 



St + fs{s,c)sx + fc{s,c)cx = 0, 
{s + a'{c))ct + fis,c))c^ = 0. 

Let U denote the state vector U = (s, c) and introduce the upper triangular matrix 

/ /. fc 



AiU) = 



\ 







/ 



s + a'(c) 



and the system Q can be read in matrix form as 

Ut + A{U) = 0. 
/ 



The eigenvalues of A are A" = fs and A"^ 



s + a 



-, with corresponding eigenvectors = (1, 0), 



(/c, A^ — A'') if < s < 1 and — (0, 1) if s = 0, 1. The eigenvalue A" may change sign whereas the 
eigenvalue A'^ is always positive. One can observe that for each c G / there exists a unique s* — s*{c) E 
(0,1) such that 

A'=(s*,c) = A"(s*,c) 

(see FiglU. For this couple (s*,c), A"^ — X'^, hence eigenvectors are not linearly independent and the 
problem is nonstrictly hyperbolic. 
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Any weak solution of (|4]l has to satisfy the Rankine-Hugoniot jump conditions given by 

f{sR,CR) - f{sL,CL) = 0-(s_R - Sl), 

CRf{sR, cr) - CLf{sL, cl) = cr{sRCR + a{cji) - slCl - a{cL)), 

where (s^, cl), {sr, cr) denote the left and right values of the couple (s, c) at a certain point of disconti- 
nuity. 

When Cfi — cl, the second equation reduces to the first equation and the speed of the discontinuity 
a is given by the first equation only. Now we are interested in the case cr, ^ c^. By combining the two 
equations ^ we may write 

(cr - CL)fisL, Cl) = a{cR - cl)sl + <j(,a{cR) - a{cL)) 

where 

f{sL,CL) ^ ) ^ — - if C^CL, 



-, aL(c) = <^ c-CL 

a'{c) if c — Cl- 



(9) 



SL +Ol(ch) ' 

Plugging this into first equation of dH), we obtain 

f^{sR + aL{cR)) = a{sL + aL{cR)) + f{sR,CR) - f{sL,CL) = f{sR,CR). 
Hence when cl^ cr the Rankine-Hugoniot condition ^ reduces to 

fisR,CR) _ f{sL,CL) _ 

SR + aL[cR) sl + cll{cr) 

3 Riemann problem 

In this section we solve the Riemann problems associated with our system, that we solve system (HI with 
the initial condition 

s(x 0) = [ < °' c(x 0) = / ^ < (10) 

^ ' ^ \ SR if a;>0 ' ^ ' ' \ CR if 2:>0 ' ^ ' 

Solution to ([Tol l is constructed by using elementary waves associated with the system. There are two 
families of waves, refered to as the s and c families, s waves consist of rarefaction and shocks (or contact 
discontinuity) across which s changes continuously and discontinuously respectively, but across which c 
remains constant, c waves consist solely of contact discontinuities, across which both s and c changes such 

that remains constant in the sense of (|9|. 

s + a'{c) 

First define a function gl by 

f a(c) - a(cL) 

- f\ ) ir c^cl, 

olIc) = < c-CL 

I a'(c) if c = Cl- 

We will restrict to the case cl > cp. The case cl < cji can be treated similarly. 

When Cl > cr the flux functions for the first equation ^ s ^ f{s, cl) and s f{s, cr) are as 
represented in Fig. |2] that is f{s,CL) < f{s,CR) Vs G (0, 1). Let Ol and 6^ be the points at which 
f{.,CL) and /{-jCr) reach their maxima respectively. 

f ( \ 

Let s* G (0, 1) be a point at which fs(s*, cl) = — r. Now draw a line through the points 

3* + aL[cR) 

(— aL(c/;), 0) and (s*, /(s*, cl) which intersects the curve /(s, cjj) at a point A> s* (see Fig. |2]l. 
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Figure 2: Two flux functions /(s, c^) and /(s, c^) with > c^. 

Our study of Riemann problems separates into two cases < s* and sl > s* which themselves 
separate into several subcases. 

• Case 1: sl < s*. 

Draw a line passing through the points {sl, /{sl, cl)) and (— ai(ci^), 0). This line intersects the 
curve /(s, cr) at points s and B (see Fig. [3]). Now we divide this into two subcases. 

• Case la: sr < B 

(a) Connect (sl, cl) to (s, c_r) by c-wave with a speed 

^ ^ f{SL,CL) ^ f{s, cr) 

Sl + CLlIcr) S + QLicR)' 

(b) Next connect (s, cr) to (s^, c^) by a s-wave, along the curve /(s, Cfl) (see Fig. O. 

For example if sj^ > s and /(s, c^) and /(s, c^) are concave functions then the solution of the 
Riemann problem is given by 

{(sl, Cl) if a; < act, 
{s,cji) if (Tct < X < ast, (11) 
{sr,Cr) if X>CFst, 

where _ _ 

f{sL,CL) f{s,CR) f{s,CR) - f(sR,Cji) 

SL + aL{cR) s + aL(cR) s - sr 

Note that < CTc < Cs. 




-ulIcr) s Sl s* a srB {sl,cl) {sr,cr) 



Figure 3: Solution of Riemann problem ( fTot with sl < s* and s/? < S. 
• Case lb: sr > B. 

Draw a line passing through the points (s/j, /(s_r, cr)) and (— aL(cij), 0). This line intersects the 
curve /(s, Cl) at a point s (see Fig. |4|. 
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(a) Connect (sl, cl) to (s, c^) by a s-wave along the curve /(s, cl). 

(b) Next connect (s, cl) to {sr, cji) by a c-wave with a speed 

f{sR,CR) f{s,CL) 



SR + UlIcr) S + ULicR)' 

For example if f{s, cl) and /(s, c^) are concave functions then the solution is given by 

{sl.cl) if X < ast, 
{s{x,t),c{x,t)) ^ { {s,cl) if (Jst < X < act 
{sr, Cr) if a; > act 



(12) 



where 



f{sR, Cr) /(s, Cl) f{s, Cl) " f{sL,CL) 

an = -, — r = -, — T: as — . 

S - SL 



SR + aL{cB.) s + aL{cR) 

Note that as < ac and {sl, cl) is connected to (s, cl) by a s-shock wave and (s, cl) is connected 
to {sr, Cfs) by a c-shock wave. 



f{s,CR) 




f{s,CL) 



X = ad 



X = act 



{s,cl) 



{sl,cl) 



{sr,cr) 



-aL{cR) SL s* A B s Sr {sl-,cl) {sr,cr) 

Figure 4: Solution of Riemann problem ( fTOb with sl < s* and sr > B. 



Case 2: sl> s*- 



Case 2a: Sfi< A . 

(a) Connect (sl, cl) to {s* ,cl) by a s-wave along the curve /(s, c^). 

(b) Connect (s*, cl) to (s, cr) by a c-wave. 

(c) Connect (s, c^) to (s^, cr) by a s-wave along the curve /(s, c^) (see Fig. |5]l. 
For example if sr < s, then the solution is given by 



where 



{s{x,t),c{x,t)) = < 



(sl,cl) 

((/.)-Hf,CL),Ci) 
((/s)-Hf,Cfl),CH) 



= /s(sL, Ci), fT2 = fsis*,CL) 



f{s*.CL) 



■cll{cr)' 



if 
if 
if 
if 
if 



X < ait, 
ait < X < a2t, 
a2t < X < a^t, 
a^t < X < ait, 
X > a^t. 



"73 = /s(s,C_r), a4 = fsisR,CR). 



Here {sl, cl) is connected to [s* , cl) by a s-rarefaction wave, (s*, cl) is connected to (s, Cij) by a 
c-shock wave and (see Fig. |5]l. If s^j > s then (s, cr) would be connected to {sr, cr) by a s-chock 
wave. 
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-oi,(ci{) Sfl s s* 6'l A Si isL,CL) (s_R,Cfl:) 

Figure 5: Solution of Riemann problem ( fTOb with > s* and s/j < A. 

Case 2b: > A 

Draw a line passing through the points {sr, /{sr, cr)) and {—aL{cR),0). This line intersects the 
curve f{s, cl) at a point s (see Fig. |6l). 




X — ad 



(s,cl) 

{sl,Cl) \ (Sfl,Cfl) 



-ahicR) s* SL A s sr {sl-,cl) {sr,cr) 

Figure 6: Solution of Riemann problem ( fTOb with < s* and s/j > ^. 



(a) Connect (s^, cl) to (s, cl) by a s-wave along the curve /(s, cl), 

(b) Next connect (s, cl) to {sr, cr) by a c-wave with a speed 

_ f{sR,CR) _ f{s,CL) 



sr + ulIcr) s + aL{cR)' 
For example if < s then the solution is given by 

{sl.cl) if X < ast, 

{s{x,t),c{x,t)) ^ ■{ {s,cl) if crst<x<act, 

{sr, Cr) if X > act, 



where 



Ksr^Cr) fis,CL) 



Hs^Cl) - f{sL,CL) 



S - SL 



(13) 



SR + aL{cR) s + aL{cR) 

Note that CTs < and (s^, cl) is connected to (s, c^) by a s-shock wave and (s, c^) is connected 
to {sr, Cr) by a c-shock wave. 
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4 Conservative finite volume scliemes for the system of polymer flood- 
ing 

Let h > and define the space grid points — ih,i £ Z and for At > define the time discretization 

points tn ~ riAt for all non-negative integer n. Let A = A numerical scheme which is in conservative 
form for equation (|4|i is given by 



i 

,n+l „"+l 



(14) 



where the numerical flux F^'-^^/2 ^l+i/2 associated with the flux functions /(s, c) and g{s, c) = 
c/(s, c), and are functions of the left and right values of the saturation s and the concentration c at 

The choice of the functions F and G determines the numerical scheme. We first present the new flux 
that we call DFLU, which is constructed as presented in the introduction. We compare it with the exact 
Riemann solver and show estimates for the associate scheme. Then we recall three other schemes to 
which to compare: the upstream mobility flux and two centered schemes, Lax-Friedrichs's and FORCE. 

4.1 The DFLU numerical flux 

The DFLU flux is an extension of the Godunov scheme that we proposed and analyze in ||2] for scalar 
conservations laws with a flux function discontinuous in space. As the second eigenvalue A"^ of the system 
is always non-negative we define 



c" FT 



(15) 



■-^1+1/2 ~ S ^'i+1/2- 

. Now the choice of the numerical scheme depends on the choice of F^^-^^^- To do so we treat c{x, t) in 
/(s, c) as a known function which may be discontinuous at the space discretization points. Therefore on 
the border of each rectangle {xi-1/2, x {tn, tn+i), we consider the conservation law: 



st+/(s,cnx = 
with initial condition s{x, 0) — s° for Xi-1/2 < x < Xi^i/2- (see Fig l4.1l l. 

t = tn+l 



(16) 







0st + /(s,c^+i), = 












t — tfl 







3^1-1/2 Xi+1/2 Xi+3/2 

Figure 7: The flux functions f{-,c) is discontinuous in c at the discretization points. 
Extending the idea of [2|,we define the DFLU flux as 

pn pDFLU ( „n n n n \ 

i+l/2 ^ yi T ^i+ItH+I) 

= min{/(min{s^ 9^}, c»), /(max{sf+i, ^^VJ, cf+i)}, 
where 6'f = argmax f{-,cf). 
Remarks: 

1) Suppose c" — Co, a constant for all z,then it is easy to see that c"^^ = cq for all i. 



(17) 
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2) Suppose s f{s, c) is an increasing function (case without gravity) then 9^ — 1 for all i and from (flTl l 
we have ^^+1/2 ^ fi^?' "-D '^he finite difference scheme ( fT4] l becomes 



= 4'-A(/(s?,cf)-/(5^i,c^i)) 



which is nothing but the standard upwind scheme. 

4.2 Comparison of the DFLU flux with the flux given by an exact Riemann solver 

Now we would like to compare the exact Godunov flux F^^^^ with our DFLU flux F^^jl^ defined by 
([TtI i. For sake of brevity we considered only the case c" > c"_^i. The opposite case can be considered 
similarly. We discuss the cases considered in section |3] 

Case la: < s* , s,+i < B. See Fig. S In this case F^,/2 = /(s,, c,) = F^/^f^ . 
Case lb: s, < s*, s,+i > B. See Fig. g] 

Then Fl, = [ {J^' ''''\ < [! where a, - /(^^ ^0 - /G^^, c,) ^ ^^^^^ ^^^^ ^ 

1 ^{s^,Ci) if CT, > S-S, 

flux gives Fj^jl^ — min{/(si, q), /(max{si+i, Ci+i)}. Therefore in this case the Godunov flux 

may not be same as the DFLU flux. 

Case 2a: s, > s*, s,+i < A. See Figg] Then 
Case 2b:si > s*, s,;+i > A. See Fig|6] 



Then F^, /o = s ' A -j, ^ where a 

*+l/2 I f{Si.Ci) if tJs > ^ S-Si 

The DFLU flux is F^^jl^ — min{/(min{si, 6'^}, q), /(max{si+i, Ci+i)}. In this case these two 

fluxes are not equal, for example when as < 0. 

One can actually observe that the Godunov flux can actually be calculated with the following compact 
formula: 
Case \: Si < s*. 

f{st,c,) if/s(s,+i,c,+i) > Oor — — > -—— -, 

Si+i + aL[Ci+i) Si + aL(Ci+i) 

min(/(si, Ci),f{si, a)) otherwise, 

, _ . . , /(Si+liCi+l) f{Si,Ci) 

where Si is given by — ^ r- = — — ^— r-. 

Si+i + aL[Ci+i) Si + aL[Ci+i) 

Case 2: Si > s*. 

/(min(s,,,6'i),c,) if /^(si+i, c,+i) > or — — > 



-^1+1/2 - Si+1 + aL{Ci+i) + aL{Ci+i} 

inm{f{si,Ci),f{s^,Ci)) otherwise , 



where Si is given by 



Si+i + aL{ci+i) s.i + aL{ci+i) 

4.3 L~ and TVD bounds for the DFLU scheme 

fig, c) 

We show first L°° bounds, and TVD bounds will follow immediately. Let M — sup{/s(s, c), Vtt}- 

s,c s + a (c) 
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Lemma 1 Let sq and cq G L°°(M, [0, 1]) be the initial data and let {s"} and {c"} be the corresponding 
solution calculated by the finite volume scheme Iil4\l using the DFLUflux ( 1751 ), ( I-/7I ). When AM < 1 then 



< < 1 for allien, 
||c"||oo < ||c"-i||oo where ||c"|U = sup, 



,n — 111 1 1 „n 1 1 l„n| 

(19) 



Proof: Since < so < 1 and hence for alH, < s° < 1. By induction, assume that (fT9] l holds for all n. 
Let 

By (fT7b.it is easy to check that if \M < 1, then H — H{si, S2, S3, ci, C2, C3) is an increasing function in 
si, S2, S3 and by the hypothesis on f,H{Q, 0, 0, ci, C2, C3) = 0, i/(l, 1, 1, ci, C2, C3) = 1- Therfore 

= i7(0,0,0,c?_i,cf,cr+i) 
< i7(l,l,l,c?_i,c?,c?+i) = l. 

This proves < s"+^ < 1. 

To prove the boundness of c, consider 

+ aic^') - c^sr - a{c-)) + A(GIVi/2 - GI^i/2) = 0. 

By adding and subtracting the term c"s"+^ to the second equation of ( fT4l i and by substituting first equation 
we can rewrite the second equation as 

where a(c^+^) - a(cf ) = a'(C^^^^)(cr^^ - ) for some ^^+^ between c^'+^ and . This is equivalent 
to 

j — 1/2 

which is the scheme written in the non-conservative form. Let 6" = A — -777; — then 

= (1 - + &rc^i < max{cr, c^l} if &r < i- 

This proves the second inequality. I 

Since c"+^ is a convex combination of c" and c"_]^ if AM < 1, then we obtain the following total 
variation diminishing property for : 

Lemma 2 Lef {c"} /je f/ze solution calculated by the finite volume scheme ( 1741 ), ( 1751 ), ( 1771 ). When XM < 1 
f/zen 

Note that the saturation itself is not TVD because of the discontinuity of /, and that the above proof applies 
also to the usptream mobility flux presented below. 

4.4 The upstream mobility flux 

Petroleum engineers have designed, from physical considerations, another numerical flux called the up- 
stream mobility flux. It is an ad-hoc flux for two-phase flow in porous media which corresponds to an 
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approximate solution to the Riemann problem. For this flux G"^^-^^^ is given again by (fTsT i and F'^^i is 
given by 



( A,(s?, cf) if (7 + - 9^)\} >0, i = l,2,t^£, 
[ A,(sr+i,cf+i) if q + {g, - g,)X*, <0, 1,2,1 j^l, 

4.5 The Lax-Friedrichs flux 

In this case fluxes are given by 

/ n _ n\ 



4.6 The FORCE flux 



This flux 11221 [3l. introduced by E. F. Toro, is an average of the Lax-Friedrichs and Lax-Wendroff flux. It is 
defined by 

Fr+i/2 = i[/(^r+i,cr+i) + /(.r,cn + 2/(.r^/^)- ^ 'i 

G'r+i/2 = i[cr+i/(sr+i,cr+i)+cr/k",c?) + 2cri/V(5r'/',cr'/') 
(c;^i«m + ««Vi)-<'<'-«(cr)), 



where 



and 



A 



n+1/2 n+1/2 /n+l/2x _ (s"+iC"_^i + s"c") ^ (a(cf+i) + a(cf )) 



i(cIVi/(sm>c?n)-<7(sr,cr)). 



5 Numerical experiments 

To evaluate the performance of the DFLU scheme we first compare its results to an exact solution and 
evaluate convergence rates, and then compare it with other standard numerical schemes already mentioned 
in the previous section, that are the Godunov, upstream mobility, Lax-Friedrichs and FORCE schemes. 

5.1 Comparison with an exact solution 

In this section we compare the calculated and exact solutions of two Riemann problems. We consider the 
following functions 

/(s,c) = s(4-s)/(l + c), a{c)=c. (20) 

Note that /(O, c) = /(4, c) = for all c and the interval for s is [0, 4] instead of [0, 1]. This choice of 
/, which does not correspond to any physical reality, was done in order to try to have a large difference 
between the Godunov and the DFLU flux (see second experiment below). 
In a first experiment the initial condition is 

I n\ / 2-5 if a; < .5, , , / .5 if x < .5, 

1 if a;>.5 ' if x>.5. ^21) 
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These / and initial data correspond to the case 2a in sections[3]and l4.2l where the DFLU flux coincides with 
the Godunov flux: F^^^^ [sl,sr,cl,cr) ^ {sl,sr,cl,cr) With s* = 1.236, A = 2.587, s = .394. 
The exact solution of the Riemann problem at a time t is given by 



s{x, t) 



2.5 if x < .5 + (71 <, 

i(4-1.5(^)) if .5 + (Tit < a; < .5 + crct 

s = .394 if .5 + (Tc i < a; < -5 + (72 t 

1. if a; > (72t + .5 



where ai ^ fs{sL,CL) = -2/3, = fs{s*,CL) 
f{s, cr) - f{sR, Cr) 



f{s*,CL) f{s,CR) 



.5 if X < .5 + (Jet, 
0. if x > .5 + act. 



s*+aL(cfl) s + aLicR) 



(22) 

1.018 and (72 = 



= 2.606. 



s -sr 



Figs, m and |9] verify that the DFLU and Godunov schemes give coinciding results. As expected both 
schemes are diffusive at c-shocks as well as at s-shocks but as the mesh size goes to zero calculated 
solutions are getting closer to the exact solution (see Fig|9]l. Table [T] shows Li errors for s and c and the 
convergence rate a. Calculations are done with ^ — j [M = 4), that is the largest time step allowed by the 
CFL condition. 
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Figure 8: Comparison with exact solution of Riemann problem 

for/i= 1/100, A = 1/4. 
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, (EB: s (left) and c (right) at i = .5 
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Figure 9: Comparison with exact solution of Riemann problem ( l20b . ( |2TI ): s (left) and c (right) at t = .5 

for/i= 1/800, A ^ 1/4. 
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h 


Godunov,||s — s/iH^i 


a 


DFLU,||s-s,,||ii 


a 


1/50 


.2373 




.2372 




1/100 


0.15134 


0.6489 


0.1506 


0.655 


1/200 


9.6868 xlO-^ 


0.6437 


9.6868 xlQ-^ 


0.6366 


1/400 


6.4228 xlQ-^ 


0.5928 


6.4228 xlQ-^ 


0.5928 


1/800 


4.2198 xlQ-^ 


0.606 


4.2197 xlQ-^ 


0.606 



h 


Goduno V, 1 1 c — c/,, | ^ i 


a 


DFLU,||c-c,,||ii 


a 


1/50 


6.3796 xlQ-^ 




6.3796 xlO-^ 




1/100 


4.1630 xlQ-^ 


0.6158 


4.1630 xlO"-' 


0.6158 


1/200 


2.6669 xlQ-^ 


0.6424 


2.6669 xlO-^ 


0.6424 


1/400 


1.7398 xlQ-^ 


0.6162 


1.7398 xlQ-^ 


0.6162 


1/800 


1.1522 xlQ-^ 


0.5945 


1.1522 xlQ-^ 


0.5945 



Table 1: Riemann problem ( |20| |. (l2Tl i: L^-errors between exact and calculated solutions at i = .5 



Now we want to have an experiment where the DFU flux differs from the Godunov flux. Therefore we 
now consider the Riemann problem with initial data 



(a;,0) 



2.3 
3.2 



if 
if 



X < 
X > 



c(j:,0) 







if 
if 



X < .5, 
X > .5. 



(23) 



This initial data corresponds to case 2b of sections [3] and l4.2l with cj^ = 0, s* = 1.236. In this case, the 
exact solution of the Riemann problem at a time t is given by 



s{x, t) 



SL = 2.3 
s = 2.7536 
SR = 3.2 



if 
if 
if 



f{sL-CL) - I{s,Cl) 



X < + CTgi 

.b + Ggt < X < .1 

X > Get + .5 

.702, and CTc 



ant. 



c(a;,0) 



if 
if 



X < 
X > 



f{sR, Cr) 



0.609. 



where (Ts _ , , 

sl-s SR + aL(CR) 

Figs. [To] and [TT] show the comparison of the results obtained with the DFU and Godunov fluxes with 
the exact solution. The solution obtained with the DFU and Godunov flux are very close even if they do 
not coincide actually. Table |2] shows Li errors for s and c and the convergence rate a. Calculations are 
done with X — j [M — 4), that is the largest time step allowed by the CFL condition. 
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Figure 10: Comparison with exact solution of Riemann problem ( l20t . ( l23l i: s (left) and c (right) at i = .5 
for/i= 1/100, A = 1/4. 
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Figure 11: Comparison with exact solution of Riemann problem ( |20] |. (l23T i: s (left) and c (right) att — .5 

foTh= 1/800, A = 1/4. 



h 


Godunov, |s— s/iH^i 


a 


DFLU,||s-s^||ii 


a 


1/50 


0.10246 




0.10373 




1/100 


5.7861 xlO-^ 


0.8243 


5.8731 xlQ-^ 


0.8206 


1/200 


3.2849 xlO-^ 


0.81674 


3.3259 xlQ-^ 


0.8203 


1/400 


1.9152 xlO-^ 


0.7785 


1.9353 xlO--' 


0.7811 


1/800 


1.1489 xlO-^ 


0.7370 


1.1571 xlQ-^ 


0.7420 



h 


Godunov, 1 1 c — c;,, | ^ i 


a 


DFLU,||c-c,,||ii 


a 


1/50 


4.8407 xlO-^ 




4.8486 xlQ-^ 




1/100 


3.0161 xlO-^ 


0.6825 


3.0201 xlQ-^ 


0.6829 


1/200 


1.9307 xlO-^ 


0.6435 


1.9328x10-^ 


0.6439 


1/400 


1.2618 xlO-^ 


0.6136 


1.2628 xlQ-^ 


0.6140 


1/800 


8.4125x10-^ 


0.5848 


8.4173 xlQ-^ 


0.5851 



Table 2: Riemann problem (|20] |. ( l23T i: L^-errors between exact and calculated solutions at t = .5. 



5.2 Comparison of the DFU, upstream mobility, FORCE and Lax-Friedrichs fluxes 

In the previous section, we have seen that Godunov and DFLU fluxes give schemes with very close per- 
formances. In this section we compare the DFLU flux with the other fluxes that we mentioned in section|4] 
which are the upstream mobility, FORCE and Lax-Friedrichs fluxes. We take now 

/(s,c) = ipi^ ^[V+i9l -32)A2(s,c)], 

Ai(s,c) + \2(S,C) 

(24) 

Ai(s,c) = — — , A2(s,c) = (1 - sY, gi = 2,gi = 1,(^ = 0, 

.5 + c 
a(c) = .25c. 

In all following experiments the discretization is such that At = 1/125 and h — 1/100. 

We first consider a pure initial value problem. Initial condition (see top of Fig. fTST i is given by 

. . / .9 if X < .5, , . / .9 if x < .5, 

s{x,\)) — < , .„ ^ , cx,0) = < „ .„ ^ . (25) 

^ ' ' [ .1 if X > .5 ' ^ ' ' \ .3 if x > .5 

With this initial condition we have F^^^^ {sl, sr, cl,cr) — F'-'{sl, sr, cr, cl) with sl — -9, s^, — 
.1,cl — 1. and cr = 0. Boundary data are such that 

s(0,t) = .9, s(2,i) = .1, c(0,i) = .9, c(2,i) = .3 V t > 0. (26) 
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Calculated solutions at time levels t=l and 1 .5 are shown in FiglT2] They show that, as expected, the DFLU 
flux, which is the closest to a Godunov scheme, performs better than the other schemes. The upstream 
mobility flux, which is an upwind scheme, performs better than the two central difference schemes, the 
FORCE and Lax-Friedrichs schemes. 
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Figure 12: s (left) and c (right) calculated at t=0., t=l. and t=1.5 for data (|24li, (|25]l, (|26] |. 

To confirm these first observations we consider now a boundary value problem. We just changed 
the boundary functions, so instead of boundary conditions dZST l we consider now a problem with closed 
boundaries, that is fluxes are zero at the boundary: 

/ = at a; = and x = 2 for ah i > 0. (27) 
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They show that, as expected, the DFLU scheme, which is the closest to a Godunov scheme, performs better 
than the upstream mobility, the FORCE or the Lax-Friedrichs schemes. 





Figure 13: s (left) and c (right) calculated at t=l., t=2. and t=3. for data dJUl, ( |26] |. 

The purpose of the last experiment whose results are shown in Fig. [T4]is to show the effect of polymer 
flooding. In this experiment we remove polymer flooding and take c = at all time. By comparing with 
the solution shown in Fig. [T3]bottom left we observe that as expected the saturation front is moving faster 
since there is no retardation due to the increase of viscosity of the wetting fluid caused by the polymer 
injection. We also observe that the structure of the solution is less complex. 
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Figure 14: s (left) calculated at t=l. and t=3. for same data as in Fig. [T3]but without polymer injection. 

6 Conclusion 

The DFLU flux defined in |2| for scalar conservation laws was used to construct a new scheme for a class 
of system of conservation laws. It was applied to a system for polymer flooding. It is very close to the 
flux given by an exact Riemann solver and the corresponding finite volume scheme compares favorably to 
other schemes using the uptream mobility, the Lax-Friedrichs and the FORCE fluxes. The DFLU is also 
very easy to implement. The extension to the case with a change of rock type is straightforward since 
the DFLU flux was built to solve this case. It will work even in cases where the upstream mobility fails 
|[T6l. In a separate paper 1 1 1 we show how to use the DFLU flux to solve Hamilton-Jacobi equations with 
a discontinuous Hamiltonian. 
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